# Figures 10 and 11

library(RColorBrewer)
library(cshapes)
#set the working directory here
setwd("")

# load the IGO membership array:
load("membership-expanded3.rda")

# load the communities data:
load("clubs of clubs data 2.rda")

# load the IGO statistics:
load("igostats5.rda")

source("ccode2cow.R")
source("report.progress.R")


# make an IGO-year dataset:
IGOlist<-dimnames(membership)$IGO
yearlist<-sort(unique(comm$year))
igos<-data.frame(IGO=rep(IGOlist, length(yearlist)), year=rep(yearlist, each=length(IGOlist)), Africa=NA, Europe=NA, LatAm=NA, Other=NA)

for (i in 1:nrow(igos)){
	report.progress(i, nrow(igos))
	igo.current<-igos$IGO[i]
	year.current<-igos$year[i]
	
	# get the full list of member states for the current IGO:
	members.current<-names(which(membership[,as.character(igo.current), as.character(year.current)]==1))
	
	# get the list of countries in each community for the current year:
	Africa.current<-comm$cow[comm$year==year.current & comm$comm.name=="Africa"]
	Europe.current<-comm$cow[comm$year==year.current & comm$comm.name=="Europe"]
	LatAm.current<-comm$cow[comm$year==year.current & comm$comm.name=="Latin America"]
	Other.current<-comm$cow[comm$year==year.current & comm$comm.name=="Other"]
	
	# count up the matches:
	igos$Africa[i]<-sum(Africa.current %in% members.current)
	igos$Europe[i]<-sum(Europe.current %in% members.current)
	igos$LatAm[i]<-sum(LatAm.current %in% members.current)
	igos$Other[i]<-sum(Other.current %in% members.current)

	}

# Now calculate the relative representation of each community in each IGO:
igos$total<-igos$Africa + igos$Europe + igos$LatAm + igos$Other
igos$Africa.prop<-igos$Africa/igos$total
igos$Europe.prop<-igos$Europe/igos$total
igos$LatAm.prop<-igos$LatAm/igos$total
igos$Other.prop<-igos$Other/igos$total

# Now calculate Herfindahl indices:
igos$Herfindahl<-igos$Africa.prop^2 + igos$Europe.prop^2 + igos$LatAm.prop^2 + igos$Other.prop^2



# Plot changes over time in each category:
threshold<-2/3 # choose a threshold for the Herfindahls (or for the proportions)
yearlist<-sort(unique(comm$year))
counts<-matrix(NA, nrow=length(yearlist), ncol=8, dimnames=list(year=yearlist, group=c("Total", "Regional", "Nonregional1", "Africa", "Europe", "LatAm", "Other", "Nonregional2")))
for (y in yearlist){
	counts[as.character(y), "Total"]<-sum(!is.na(igos$Herfindahl) & igos$year==y, na.rm=T)
	counts[as.character(y), "Regional"]<-sum(igos$Herfindahl>=threshold & igos$year==y, na.rm=T)
	counts[as.character(y), "Nonregional1"]<-sum(igos$Herfindahl<threshold & igos$year==y, na.rm=T)
	counts[as.character(y), "Africa"]<-sum(igos$Africa.prop>=threshold & igos$year==y, na.rm=T)
	counts[as.character(y), "Europe"]<-sum(igos$Europe.prop>=threshold & igos$year==y, na.rm=T)
	counts[as.character(y), "LatAm"]<-sum(igos$LatAm.prop>=threshold & igos$year==y, na.rm=T)
	counts[as.character(y), "Other"]<-sum(igos$Other.prop>=threshold & igos$year==y, na.rm=T)
	counts[as.character(y), "Nonregional2"]<-counts[as.character(y), "Total"]- sum(counts[as.character(y), c("Africa", "Europe", "LatAm", "Other")])
	}
	
colorscheme<-brewer.pal(5, "Spectral")
names(colorscheme)<-c("Other", "Latin America", "Non-regional", "Europe", "Africa") 

pdf("Figure 10.pdf", width=7, height=7)
plot(yearlist, yearlist, type="n", las=1, bty="n", xlab="Year", ylab="Number of IGOs", ylim=c(0,400))
polygon(c(yearlist, rev(yearlist)), c(counts[,"Africa"], rep(0, length(yearlist))), lwd=0.1, col=colorscheme["Africa"])
polygon(c(yearlist, rev(yearlist)), c(counts[,"Africa"], rev(counts[,"Africa"]+counts[,"LatAm"])), lwd=0.1, col=colorscheme["Latin America"])
polygon(c(yearlist, rev(yearlist)), c(counts[,"Africa"]+counts[,"LatAm"], rev(counts[,"Africa"]+counts[,"LatAm"]+counts[,"Other"])), lwd=0.1, col=colorscheme["Other"])
polygon(c(yearlist, rev(yearlist)), c(counts[,"Africa"]+counts[,"LatAm"]+counts[,"Other"], rev(counts[,"Africa"]+counts[,"LatAm"]+counts[,"Other"]+counts[,"Europe"])), lwd=0.1, col=colorscheme["Europe"])
polygon(c(yearlist, rev(yearlist)), c(counts[,"Africa"]+counts[,"LatAm"]+counts[,"Other"]+counts[,"Europe"], rev(counts[,"Africa"]+counts[,"LatAm"]+counts[,"Other"]+counts[,"Europe"]+counts[,"Nonregional2"])), lwd=0.1, col=colorscheme["Non-regional"])
text(2004, 50, "African", adj=1)
text(2004, 110, "Latin American", adj=1)
text(2004, 144, "South Asian/Other", adj=1)
text(2004, 190, "European/Northern", adj=1)
text(2004, 280, "Non-cluster", adj=1)
dev.off()


### Now identify the most central countries in each community:

# Refer to Tore Opsahl's tnet package for dealing with weighted networks:  http://toreopsahl.com/tnet/weighted-networks/node-centrality/
library(tnet)

# first make an array containing the number of shared IGO memberships across dyads (this matches what was used to generate the communities in the first place)
shared<-array(NA, dim=c(dim(membership)[1], dim(membership)[1], dim(membership)[3]), dimnames=list(dimnames(membership)[[1]], dimnames(membership)[[1]], dimnames(membership)[[3]]))
for (i in 1:dim(shared)[3]){
	shared[,,i]<-membership[,,i] %*% t(membership[,,i])
	diag(shared[,,i])<-NA
	}

year.current<-2005
commlist<-sort(unique(comm$comm.name))
mostcentral<-rep(NA, length(commlist))
names(mostcentral)<-commlist
for (i in 1:length(commlist)){
	comm.current<-commlist[i]
	memberlist<-comm$cow[comm$year==year.current & comm$comm.name==comm.current]
	adj.current<-shared[memberlist, memberlist, as.character(year.current)]

	# Now convert it into the edgelist format for use by the tnet package:
	source("matrix2dyads.R")
	edgelist.current<-matrix2dyads(adj.current)
	edgelist.current<-edgelist.current[edgelist.current[,1]!=edgelist.current[,2],] # remove self-loops
	
	# Create sequential node numbers for use by the tnet package (note that I can't use the COW IDs here because they're not sequential)
	cowlist<-sort(unique(c(edgelist.current[,1], edgelist.current[,2])))
	nodelist<-1:length(cowlist)
	edgelist.current[,1]<-nodelist[match(edgelist.current[,1], cowlist)]
	edgelist.current[,2]<-nodelist[match(edgelist.current[,2], cowlist)]

	# Now do the degree centrality calculation and replace the node numbers with the COW labels:
	net.current<-as.tnet(edgelist.current, type="weighted one-mode tnet")
	a<-degree_w(net.current)
	rownames(a)<-cowlist[match(a[,1], nodelist)]
	mostcentral[comm.current]<-rownames(a)[which.max(a[,"output"])]
	}

mostcentral # SEN, FRN, VEN and IND for 2005



### Figure 11:

threshold<-2/3 # choose a threshold for the proportions
countrynames<-read.csv("COW country codes.csv")
pdf(file="Figure 11.pdf", width=8, height=8)
par(mfrow=c(2,2))
for (cow.current in mostcentral){
	m<-matrix(NA, nrow=length(yearlist), ncol=6, dim=list(year=yearlist, group=c("Total", "Africa", "Europe", "LatAm", "Other", "Nonregional")))
	m[,"Total"]<-colSums(membership[cow.current,,as.character(yearlist)])
	for (y in yearlist){
		IGOs.current<-names(which(membership[cow.current,,as.character(y)]==1))
		m[as.character(y), "Africa"]<-sum(IGOs.current %in% igos$IGO[igos$Africa.prop>=threshold & igos$year==y & !is.na(igos$Africa.prop)])
		m[as.character(y), "Europe"]<-sum(IGOs.current %in% igos$IGO[igos$Europe.prop>=threshold & igos$year==y & !is.na(igos$Europe.prop)])
		m[as.character(y), "LatAm"]<-sum(IGOs.current %in% igos$IGO[igos$LatAm.prop>=threshold & igos$year==y & !is.na(igos$LatAm.prop)])
		m[as.character(y), "Other"]<-sum(IGOs.current %in% igos$IGO[igos$Other.prop>=threshold & igos$year==y & !is.na(igos$Other.prop)])
		}
	m[,"Nonregional"]<-m[,"Total"]-rowSums(m[,c("Africa","Europe","LatAm","Other")])
	plot(yearlist, yearlist, type="n", las=1, bty="n", xlab="Year", ylab="Number of IGOs", ylim=c(0,max(m)), main=countrynames$StateNme[match(cow.current, countrynames$StateAbb)])
	polygon(c(yearlist, rev(yearlist)), c(m[,"Africa"], rep(0, length(yearlist))), lwd=0.1, col=colorscheme["Africa"])
	polygon(c(yearlist, rev(yearlist)), c(m[,"Africa"], rev(m[,"Africa"]+m[,"LatAm"])), lwd=0.1, col=colorscheme["Latin America"])
	polygon(c(yearlist, rev(yearlist)), c(m[,"Africa"]+m[,"LatAm"], rev(m[,"Africa"]+m[,"LatAm"]+m[,"Other"])), lwd=0.1, col=colorscheme["Other"])
	polygon(c(yearlist, rev(yearlist)), c(m[,"Africa"]+m[,"LatAm"]+m[,"Other"], rev(m[,"Africa"]+m[,"LatAm"]+m[,"Other"]+m[,"Europe"])), lwd=0.1, col=colorscheme["Europe"])
	polygon(c(yearlist, rev(yearlist)), c(m[,"Africa"]+m[,"LatAm"]+m[,"Other"]+m[,"Europe"], rev(m[,"Africa"]+m[,"LatAm"]+m[,"Other"]+m[,"Europe"]+m[,"Nonregional"])), lwd=0.1, col=colorscheme["Non-regional"])
	if (cow.current=="SEN"){
		text(2000, 15, "African", adj=1)
		text(2000, 45, "Non-cluster", adj=1)
		text(1974, 32, "European/Northern", adj=1)
		arrows(1974.5, 32, 1980, 25, length=0.1, lwd=2)
		}
	if (cow.current=="FRN"){
		text(2000, 90, "Non-cluster", adj=1)
		text(2000, 40, "European/Northern", adj=1)
		text(1990, 4, "African", adj=1)
		}
	if (cow.current=="VEN"){
		text(2000, 10, "Latin American", adj=1)
		text(2000, 50, "Non-cluster", adj=1)
		text(1975, 25, "European/Northern", adj=1)
		arrows(1976, 25, 1982, 24, length=0.1, lwd=2)
		}	
	if (cow.current=="IND"){
		text(2000, 50, "Non-cluster", adj=1)
		text(1994, 12, "European/Northern", adj=1)
		text(1963, 4, "African", adj=1)
		text(2004, 12, "South", adj=1)
		text(2004, 7, "Asian", adj=1)
		}
	} # close country loop
dev.off()

